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Abstract 

During a crossover via a switching mechanism from one 2-body poten- 
tial to another as might be applied in modeling (chemical) reactions in 
the vicinity of bond formation, energy violations would occur due to finite 
step size which determines the trajectory of the particles relative to the 
potential interactions of the unbonded state by numerical (e.g. Verlet) 
integration. This problem is overcome by an algorithm which preserves 
the coordinates of the system for each move, but corrects for energy dis- 
crepancies by ensuring both energy and momentum conservation in the 
dynamics. The algorithm is tested for a hysteresis loop reaction model 
with an without the implementation of the algorithm. The tests involve 
checking the rate of energy flow out of the MD simulation box; in the 
equilibrium state, no net rate of flows within experimental error should 
be observed. The temperature and pressure of the box should also be 
invariant within the range of fluctuation of these quantities. It is demon- 
strated that the algorithm satisfies these criteria 
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1 PRELIMINARIES 

The dimeric particle reaction simulated may be written 

ki 

2A t± A 2 (1.1) 
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where k\ is the forward rate constant and k-% is the backward rate constant. 
The reaction simulation was conducted at extremely high temperatures which 
are off-scale and not used in ordinary simulations of L J (Lennard-Jones) fluids 
where normally pQ the reduced temperatures T* ranges ~ 0.3 — 1.2, whereas 
here, T* ~ 8.0 — f6.0, well above the supercritical regime of the LJ fluid At 
these temperatures, the normal choices for time step increments do not obtain 
without also taking into account energy-momentum conservation algorithms in 
regions where there are abrupt changes of gradient P3[2J[3]. The global literature 
does not seem to cover such extreme conditions of simulation with discrete time 
steps using the Verlet velocity algorithm. The units used here are reduced LJ 
ones pp. The simulation was at density p = 0.70 with 4096 atomic particls which 
could react. The potentials used are as given in Fig. (pp) where r& = 1.20 for the 
vicinity where the bond of the dimer is broken and where 2 free particles emerge, 
and rf = 0.85 is the point along the hysteresis potential curve where the dimer 
is defined to exist for two previously free particles which collide. The reaction 
proceeds as follows; all particles interact with the splined LJ pair potential ulj 
except for the dimeric pair formed from particles i and j which interact 
with a harmonic-like intermolecular potential modified by a switch u(r) given 

by 

u(r) = u mb (r)s(r) + u LJ [1 - s(r)] (1.2) 
where u v n,(r) is the vibrational potential given by ea. (|1.3|) below 

u V ib(r) = u + -k(r - r ) 2 (1.3) 
The switching function s(r) is defined as 

s(r) = 7~ \Tl (1.4) 

where 

s(r) ->• 1 if r < r sw 
s(r) — * for r > r sw 

The switching function becomes effective when the distance between the atoms 
approach the value r sw (see Fig. (fTJ) ). Some of the other parameters used in 
the equations that follow include: 

mo = —10, ro = 1.0, k ~ 2446 (exact value is determined by the other input pa- 
rameters), n = 100, rt = 0.85, = 1.20, and r sw = 1.11. Particles i and j above 
also interact with all other particles not bonded to it via u^j . Full simulation 
details are given elsewhere [2]; suffice to say the activation energy at ry is ex- 
tremely high at approximately 17.5. At rf, the molecular potential is turned on 
where at this point there is actually a crossing of the potential curves although 
the gradients of the molecular and free ulj potentials are "'very close"'. On 
the other hand, at r& , the switch forces the two curves to coalesce, but detailed 
examination shows that there is an energy gap of about the same magnitude as 
the cut-off point in a normal non-splincd LJ potential (~ 0.04 energy units), 
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meaning there is no crossing of the potentials. The current algorithm is applied 
for both these cross-over regions with their different mechanisms of cross-over. 
The MD cell is rectangular, with unit distance along the axis ( x direction) 
of the cell length, whereas the breadth and height was both 1/16, implying a 
thin pencil-like system where the thermostats were placed at the ends of the 
MD cell, and the energy supplied per unit time step St at both ends of the cell 
(orthogonal to the x axis) in the vicinity of x = and x = 1 maintained at 
temperatures T% and 7} could be monitored, where this energy per unit step 
time is respectively th and e/. At equilibrium, (when Th = TJ), the net energy 
supplied within statistical error (meaning 1-3 units of the standard error of the 
e distributions ) is zero, i.e. e; « eh 0. The cell is divided up uniformly into 
64 rectangular regions along the x axis and its thermodynamical properties of 
temperature and pressure are probed. The resulting values of the e's and the 
relative invariance of the pressure and temperature profiles would be a measure 
of the accuracy of the algorithm from a thermodynamical point of view at the 
steady state. For systems with a large number of particles such thermodynam- 
ical criteria are appropriate. The synthetic thermostats now frequently used 
in conjunction with " 'non-Hamiltonian" ' MD 3 cannot be employed for this 
type of study, where actual energy increments are sampled. The runs were for 
4 million time steps, with averages taken over 100 dumps, where each variable 
is sampled every 20 time steps. The final averages were over the 20-100 dump 
values of averaged quantities. 
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Figure 1: Potentials used for this work. 



The temperature T and pressure p are computed by the equipartition and 
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Figure 2: Temperature profile across the cell for different set conditions a — e for 
temperature T* and step time 5t pairs (T* ,St) where a — (8.0, 2.0 ep — 3), b = 
(8.0, 5.0 ep - 4), c = (8.0, 5.0 ep-5),d= (12.0, 5.0 ep - 5), e = (16.0, 5.0 ep - 5). 
The curves {11, 13, tl, t2, t3} results with the application of the algorithm at r& 
and r/ with associated conditions 11 a, 13 <^ b, tl c, t2 rf, t3 e whilst 
the curves {12, 14, 15, 16, 17} are for the cases without implementing the algorithm 
with the associated conditions 12 a, 14 6, 15 4=> c, 16 <^> rf, 17 e, where 
epx = 10 x . 
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Figure 3: Pressure profile across the cell for different runs. The conditions of the 
runs and the labeling of the curves are exactly as in Fig. ([2|). 
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Curve 


e/i ei Mean Temperature 


11 


-.2274E+00 ±0.19E-02 -.2295E+00 ±0.21E-02 0.9063E+01 ±0.62E-02 


12 


-.5602E+00 ±0.22E-02 -.5596E+00 ±0.22E-02 0.1032E+02 ±0.63E-02 


13 


-.4161E-01 ±0.14E-02 -.4089E-01 ±0.14E-02 0.8774E+01 ±0.79E-02 


14 


-.5201E-01 ±0.16E-02 -.5103E-01 ±0.17E-02 0.8980E+01 ±0.98E-02 


tl 


-.5312E-03 ±0.92E-03 -.3334E-03 ±0.76E-03 0.8082E+01 ±0.49E-02 


15 


0.1311E-02 ±0.82E-03 0.1147E-02 ±0.84E-03 0.7731E+01 ±0.97E-02 


t2 


-.6823E-03 ±0.12E-02 -.1507E-02 ±0.13E-02 0.1216E+02 ±0.17E-01 


16 


0.7291E-02 ±0.12E-02 0.6343E-02 ±0.14E-02 0.1088E+02 ±0.15E-01 


t3 


-.9348E-03 ±0.18E-02 -.3379E-02 ±0.17E-02 0.1622E+02 ±0.18E-01 


17 


0.1918E-01 ±0.14E-02 0.1938E-01 ±0.16E-02 0.1329E+02 ±0.20E-01 



Table 1: Table with values for the mean heat supply per unit step and temper- 
ature. The error is one unit of standard error for the quantities. 

Virial expression given respectively by 

/^Pi-Pi/mA = 3Nk B T and P = pk B T + W/V 

where W — — g ^2j>i w ( r ij) an( i the intermolecular pair Virial w(r) is given 
by w(r) = r dv ^ with v being the potential. 

2 ALGORITHM AND AND ANALYSIS OF NU- 
MERICAL RESULTS 

The velocity Verlet algorithm [4, p. 81]used here pQ and allied types generate 
a trajectory at time n5t from that at (n — l)St with step increment St through 
a mapping T m where (v(nSt),r(nSt)) — T m (v((n — l)St),r((n — l)St)) which 
does not scale linearly with St. For a Hamiltonian H whose potential V is 
dependent only on position r having momentum components pi, the system 
without external perturbation has constant energy E, and the normal assump- 
tion in MD (NAMD)is that for the n th step, AE n = \H(nSt) - E\ < e and also 
J2iLi < e s for the specified e's. In the simulation under NAMD, the force 
fields are constant and do not change for any one time step. In these cases, the 
energy is a constant of the motion for any time interval dtr when no external 
perturbations (e.g. due to thermostat interference) are impressed. When there 
is a crossing of potentials at such a time interval interval from to <p a at an 
inter particle distance(icd) r c (such as points r $ and r& of Fig. |T])) of general 
particle 1 and 2 (the (1, 2) particle pair) due to a reactive process (such as oc- 
curs in either direction of a bifurcation occurs where the MD program 
computes the next step coordinates as for the unreacted system (potential (pf,), 
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which needs to be corrected. Let the icdat time step i be Ti (with potential) 
and at step % + 1 after interval St be rj = where r/ < r c < 7*j. Due to 
this crossover, a different Hamiltonian if' is operative after point r c is crossed, 
where under NAMD, the other coordinates not undergoing crossover are not 
affected. For what follows, subscripts refer to the particle concerned. Let the 
interparticle potential at r/ be E a = Ef — 4> a {r}) and at rf be E b = 4>(,(rf), 
where A — E b — E a . Then if r f be the final coordinate due to the 4>b potential 
and force field, two questions may be asked: (i) Can the velocities of (1,2) be 
scaled, so that there is no energy or momentum violation during the crossover 
based on the (f>b trajectory calculation and (ii) Can a pseudo stochastic potential 
be imposed from coordinates r c (at virtual time t c ) to rj such that (i) above is 
true? For (ii) we have 

Theorem 2.1 A virtual potential which scales velocities to preserve momentum 
and energy can be constructed about region r c . 

proof The external work done 5W on particles 1 and 2 over the time step is 
proportional to the distance traveled since these forces are constant and so for 
each of these particles i, F ex t,i-Ari = SWi where An is the distance increment 
during at least part of the time step from r c to rj. For the non-reacting trajec- 
tory over time XSt (A < 1) (virtual because it is not the correct path due to the 
crossover at r c ), 

5W 2 + SWi - {(P b {r f ) - <j) b {r c )) = A^K.E.) (2.1) 

where AJ2(K.E.) is the change of kinetic energy for the (1,2) pair from the 
First Law between the end points r/, r c . Now over time interval t c to tf, for 
the reactive trajectory, we introduce a '"virtual potential"' V mr that will lead 
to the same positional coordinates for the pair at the end of the time step with 
different velocities than for the non-reactive transition leading to the transition 

SWi + SWi - (V mr (r f ) - V mr {r c )) = A^ '{K.E.) (2.2) 

where A ^ '(K.E.) is the change of kinetic energy for the pair with V mr turned 
on and along this trajectory, the change of potential for V mr is equated to the 
change in the K.E. of the pair as given in the results of theorem (I2.2j) for all 
three orthogonal coordinates, i.e. 

SV mr (r) - 6<p b (r) = 6 ( A ^; /v ./•.., „ . i - '(*.£.)*,„,«) 

with momentum conservation, that is (5y 4r (r i ) = 6<fi a (ri) for the variation 
along the i\ coordinate, but 54> a (ri) = —SK.E. along internuclear coordinate fj 
whereas SV mr = —K.E. (scaled about all three axes). Continuity of potential 
implies 

Mrf) = V vir (r f );Mrc) = V vir (r c y,<p b (rc) = V vir (r c y, (2.3) 
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Subtracting (|2.1[) from (|2.2p and applying b.c.'s (|2.6p leads to 

A = ^(r/) - V^fr,) = Mrf)-<l>a(rf)=E b -E a (2.4) 

A A £(*J5.) (2.5) 

The above shows that a conservative virtual potential could be said to be oper- 
ating in the vicinity of the transition (from t c tot a ) .• 
Question (i) above leads to: 

Theorem 2.2 Relative to the velocities at any r/ due to the <pb potential, the 
rescaled velocities v ' due to the potential difference A leading to these final 
velocities due to the virtual potential can have a form given by 

v; ' = (1 + a) Vi + (3 (2.6) 

(where i = 1,2) for a vector f3. 

proof From the v velocities at r/ due to <\>\, we compute the v' velocities at 
rf due to the virtual potential. Since net change of momentum is due to the 
external forces only, which is invariant for the (1, 2) pair, conservation of total 
momentum relating v' and v in (|2.6p yields a definition of (3 ( summation from 
1 to 2 for what follows, where the mass of particle i is rrij) 

/3 = -a ^ mjVi/ ruj (2.7) 

Defining for any vector s 2 = s.s,/3 2 = a 2 Q, where 

Q= (E m * vi ) 2/ (E m *) 2 ( 2 -8) 

then the rescaled velocities become from (|2.6p 

vf = (1 + a) 2 v ; 2 + 2(1 + a)vi./3 + f3 2 . (2.9) 
With A = Eb — E a , Energy conservation implies 

E^ v S 2 -E^ v = 2 = A ( 2 - 10 ) 

The coupling of (|2. 9112.10]) leads, after several steps of algebra to 

A=- Vl 2 + v 2 2 -2 Vl .v 2 2.11 

2(mi + m 2 ) 

2am 1 m 2 r 2 2 -, 
+ ^7 : r vi +v 2 -2vi.V 2 . 

2(mi + m 2 ) 

Defining a = (vi — v 2 ) 2 , q = m 2 mi/[2(mi +m 2 )], (g > 0, a > 0), then the 
above is equivalent to the quadratic equation 

a 2 qa + 2qaa - A = (2.12) 
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and in simulations, only a is unknown and can be determined from (|2.12p where 
real solutions exist for A/qa > — 1. • The above Inequality leads to a certain 
asymmetry concerning forward and backward reactions, even for reversible re- 
actions where the region of formation and breakdown of molecules are located 
in the same region with the reversal of the sign of approximate A. For this 
simulation, a reaction in either direction (formation or breakdown of dimer ) 
proceeds if (??) is true; if not then the trajectory follows the one for the initial 
trajectory without any reaction (i.e. no potential crossover). 
Interpretation of results. Fig. |T]) shows a rapidly changing potential curve 
with several inflexion points used in the simulation at very high temperature 
(as far as I know such ranges have not been reported in the literature for non- 
synthetic methods) warranting smaller time steps; larger ones would introduce 
errors due to the rapidly changing potential and high K.E.; thus, even with 
the application of the algorithm between cordinates rj and r&, curves 11 and 
12 have too large a St value to achieve equilibrium - meaning flat or invariant 
- temperature (see Fig. ) or pressure (see Fig. d3]))or unit step thermostat 
heat supply (see Table [T])(e^ and ej) profiles where for these curves, the (e/j, e;) 
values show net heat absorption; the curve at tl (with St = 5.0 ep — 5 show flat 
profiles (within statistical fluctuations and 2 standard errors of variation) for 
temperature, pressure and net zero heat supply; and this choice of time step 
interval was found adequate for runs at much higher temperatures (T = 12 and 
T = 16) which was used to determine thermodynamical properties [2]. For this 
St value and all others, no reasonable stationary equilibrium conditions could 
be obtained without the application of the algorithm (curves 12,14,15,16 and 17). 
The algorithm is seen to be effective over a wide temperature range for this 
complex dimer reaction simulated under extreme values of thermodynamical 
variables and the results here do not vary for longer runs and greater sampling 
statistics (e.g. 6 or 10 million time steps). The thin, pencil-like geometry of the 
rectangular cell with thermostats located at the ends would highlight the energy 
non-conservation leading to a non-flat temperature distribution, as observed and 
which was used to determine the regime of validity of the algorithm. 
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